//******************************************************************************
// Date: 2020-4-14
// Author: Che Yuzhi
// ID: 201928000907028
// Description:
// 	Homework3
// 	1. 画出一个直线形状的本底和两个相邻的高斯信号峰的直方图 (用TH1F)，并且
// 	画出它们叠加后的直方图。
//	2. 拟合总的直方图，在图中显示拟合后的直方图的参数。
//	3.将图片存储为eps图片，和C++脚本一起放在作业中。 (总的统计量为10000，本
//	底和信号的参数可以自己设定。要求高斯信号峰间有重叠的部分。)
//******************************************************************************
{
void setStyle();
setStyle();

/********  Sum of statistic 10000  ********/
int signal_1 = 4250;
int signal_2 = 2750;
int background = 3000;

TF1* f1 = new TF1("f1","TMath::Gaus(x,8,1.6)",0,20);
TF1* f2 = new TF1("f2","TMath::Gaus(x,12,1.6)",0,20); // First two Gaus Fuctions
TF1* f3 = new TF1("f3","x",0,20); // Linear Fuction

int nbin = 80;
TH1F* h1 = new TH1F("h1","Signal_1",nbin,0,20);
h1->FillRandom("f1",signal_1);
TH1F* h2 = new TH1F("h2","Signal_2",nbin,0,20);
h2->FillRandom("f2",signal_2);
TH1F* h3 = new TH1F("h3","Background",nbin,0,20);
h3->FillRandom("f3",background);

TH1F* hsum = (TH1F* )h1->Clone("hsum"); // To contain all other histograms
hsum->Add(h2); 
hsum->Add(h3); // Define h1 + h2 + h3


//----------------------------------------------------
// Fit histogram using Root
//----------------------------------------------------
TF1* fitter = new TF1("fitter",
			  "[0]*TMath::Gaus(x,[1],[2])+[3]*TMath::Gaus(x,[4],[5])+[6]*x",
			  0,20); // Function to fit
fitter->SetParameters(4200,8,1.5,2750,12,1.5,3000);
fitter->SetParNames("s1-Max","s1-Mean","s1-RMS",
										 "s2-Max","s2-Mean","s2-RMS",
										 "Background-k");
fitter->SetParLimits(1,6,10);
fitter->SetParLimits(2,1,2);
fitter->SetParLimits(4,10,14);
fitter->SetParLimits(5,1,2);
hsum->Fit(fitter,"0"); 
	// Option 0: "0" Do not plot the result of the fit. 
	//  By default the fitted function is drawn 
	//  unless the option"N" above is specified.
	// Setting "N" is also OK!

TCanvas* c1 = new TCanvas();
// Set style
h1->SetLineColor(kBlue);
h2->SetLineColor(kPink);
h3->SetLineColor(kBlack);
hsum->SetLineColor(kGreen);

// Draw
hsum->Draw();
fitter->Draw("same");
h1->Draw("same");
h2->Draw("Same");
h3->Draw("Same");
c1->SaveAs("Result.Root.eps");

// Fit histogram using RooFit
//----------------------------------------------------
// One can import a TH1F to a RooDataHist belonging to  RooFit.
// And use PDF in RooFit to fit it.
// Look: https://root.cern.ch/doc/master/rf102__dataimport_8C.html
//----------------------------------------------------

// Import hsum into a RooDataHist
RooRealVar x("x","x",0,20);
RooDataHist hsum_Roo("hsum_Roo","hsum_Roo",x,RooFit::Import(*hsum));
RooPlot* frame = x.frame(RooFit::Title("Fit hsum in RooFit"));
hsum_Roo.plotOn(frame);

// Create RooFit PDF to Fit
RooRealVar mean1("mean1","mean1",7,0,10);
RooRealVar sigma1("sigma1","sigma1",1.6,1,2);
RooGaussian gauss1("guass1","gauss1",x,mean1,sigma1); // First gauss signal

RooRealVar mean2("mean2","mean2",12,10,20);
RooRealVar sigma2("sigma2","sigma2",1.6,1,2);
RooGaussian gauss2("guass2","gauss2",x,mean2,sigma2); // Second gauss

RooGenericPdf bac("bac","bacground","x", RooArgSet(x)); // bacground Linear func

RooRealVar par1("par1","par1",0.5,0,1);
RooRealVar par2("par2","par2",0.5,0,1);
RooAddPdf fitter_Roo("fitter_Roo","Sum of 3 PDFs",RooArgList(gauss1,gauss2,bac),RooArgList(par1,par2));
	// Define: signal1pdf + signal2pdf + backgroundpdf = Total PDF for fitting

// Fit the imported data
fitter_Roo.fitTo(hsum_Roo);
fitter_Roo.plotOn(frame);
fitter_Roo.paramOn(frame);

//
TCanvas* c2 = new TCanvas();
frame->Draw();
h1->Draw("SameE1");
h2->Draw("SameE1");
h3->Draw("SameE1");
c1->SaveAs("Result.Root.eps");
}

void setStyle(){
	gStyle->SetOptStat(1111);
	gStyle->SetOptFit(1111);  
}

